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ABSTRACT 

In order to study the origin of the architectures of low mass planetary systems, 
we perform numerical surveys of the evolution of pairs of coplanar planets in the 
mass range (1 — 4) M®. These evolve for up to 2 x 10 7 yr under a range of orbital 
migration torques and circularization rates assumed to arise through interaction with 
a protoplanetary disc. 

Near the inner disc boundary, significant variations of viscosity, interaction with 
density waves or with the stellar magnetic field could occur and halt migration, but 
allow circularization to continue. This was modelled by modifying the migration and 
circularization rates. 

Runs terminated without an extended period of circularization in the absence of 
migration torques gave rise to either a collision, or a system close to a resonance. 
These were mostly first order with a few % terminating in second order resonances. 
Both planetary eccentricities were small <0.1 and all resonant angles liberated. This 
type of survey produced only a limited range of period ratios and cannot reproduce 
Kepler observations. 

When circularization alone operates in the final stages, divergent migration occurs 
causing period ratios to increase. Depending on its strength the whole period ratio 
range between 1 and 2 can be obtained. A few systems close to second order com- 
mensurabilities also occur. In contrast to when arising through convergent migration, 
resonant trapping does not occur and resonant angles circulate. Thus the behaviour 
of the resonant angles may indicate the form of migration that led to near resonance. 

Key words: planetary systems: formation planetary systems: protoplanetary discs 
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1 INTRODUCTION 

The Kepler space telescope has been very successful in re¬ 
vealing the diversity of the architectures of exoplanetary sys¬ 
tems. To date, Kepler has de tected over 4000 planet candi¬ 
dates Isee lBatalha et all2013ll . from which about a thousand 
are confirmed planets. About one third of Kepler’s candi¬ 
dates are associated with multiple transiting systems and 
the vast majority of them are expected to be multi-planetary 
systems dLissauer et al.|[2012l l. Amongst these, Super-Earths 
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with orbital periods less than 20 days are abundant around 
Sun-like stars. It is unlikely that these planets formed at 
their current locations, but instead they could have mi¬ 
grated iuwards_^orn_significantly larger radii (see for ex¬ 
ample [lzidoro^t((al]|201J, for a discussion). 

Due to interaction between the planets and the proto¬ 
planetary disc, torques that can lead to orbital migration 
are expected to occur as long as the planets a re embedded 
in the disc le.g. iGoldreich fe SchlichtineH2014r i. For planets 
with different masses, migration rates are expected to differ 
leading to a time-dependent evolution of their orbital pe¬ 
riod ratio. When the outermost planet is more massive it 
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migrates through the disc faster leading to convergent mi¬ 
gration. If the inner planet is more massive, as long as its 
migration stalls at the inner boundary of the disc, conver¬ 
gent migration will then ultimately ensue. 

For systems in which the planetary orbits converge, 
the planets tend to enter and become locked in a mean- 
motion resonance and migrate together su bsequently (e.g. 
iNelson fe Papaloizoull2002 : Kiev et aP[2004l 'l. In the simplest 
case of nearly circular and coplanar orbits, first-order reso¬ 
nances are readily formed between two planets undergoing 
convergent migration. For these the orbital period ratio can 
be expressed as (p + 1 )/p, with p being an integer. 

Kepler’s candidate multi-planetary systems show indi¬ 
cations of having been affected by resonances. Notably there 
are excesses in the numbers of sys tems .just wid e of the 2:1 
and 3:2 commensurabilities fe.g. Isteffenl 120131 1. However, 
most systems are non resonant with a wide range of pe¬ 
riod ratios for consecutive planets being observed. This is 
apparently in conflict with the notion of widespread conver¬ 
gent migration. However, much of the previous discussion 
of this issue has been based on utilisation of simplified pro¬ 
toplanetary disc models without consideration of conditions 
near the inner boundary where migration can be halted. 
The proper incorporation of such effects is problematic as it 
is likely to require detailed modelling of structural features 
in the protoplanetary disc resulting from, either significant 
variations of effective viscosity if the transition is between 
regions of varying levels of turbulence, or the interaction of 
the disc with the magnetic field of the central star. In the 
latter case there may be a region where the surface density 
increases outwards supporting outward propagating density 
waves. In this paper we assess the effects of a transition at an 
inner disc boundary, or the possible effects of wake-planet or 
density wave - planet interaction, on the inward migration 
of a pair of planets through simplified prescriptions for the 
orbital migration and circularization rates of the planets. 

In order to investigate the outcome of the migration of a 
pair of low mass planets and the consequent architectures of 
low mass planetary systems, we have performed numerical 
IV-body simulations of a pair of coplanar migrating low-mass 
planets in the mass range (1 — 4) M®. In the surveys we con¬ 
ducted the planets were initialised on circular coplanar or¬ 
bits with the outer one being at 1 AU or 5 AU and the inner 
one being just outside either the 2:1 commensurability or the 
3:2 commensurability. The orbital evolution was considered 
for a large range of migration and circularization times which 
might be supposed to arise from interacting with the proto¬ 
planetary disc for times up to 2 x 10 7 yr. We suppose the 
disc to possess an inner cavity where the structure changes 
on account of for example, rapid changes in disc viscosity, 
or interaction with the central star. We suppose that on ac¬ 
count of structural changes, inward migration ceases inside 
the cavity while investigating cases for which orbital circu¬ 
larization ceases for both planets, it continues for only the 
outer planet, or it continues for both planets. The last two 
cases result in divergent migration for the two planets, so 
modelling the ef fect of wake-planet, or p lanet- density wave 
interactions fsee l Baruteau fe Papaloizoull2013l . for a discus¬ 
sion). Such runs have to be considered if the final range of 
period ratios is to contain those found in the observations. 

With this modelling, we find that the runs in our sur¬ 
veys can have a variety of ends. These include a collision 


Resonance 

System 

Period ratio 

5:3 

Kepler 365 

1.66754 

5:3 

Kepler 262 

1.67322 

5:3 

Kepler 87 

1.66671 

7:5 

HD 41248 

1.394 ± 0.005 

9:7 

Kepler 29 

1.28567 

9:7 

Kepler 417 

1.29292 


Table 1. Apart from the case of HD 41248 (see text) , the en¬ 
tries are confirmed Kepler systems that are close to second-order 
commensurabilities. 


between the planets, a configuration in or close to a first 
order resonance, and systems occupying a range of period 
ratios as is found for the observed systems. In a relatively 
small number of cases, second order commensurabilities are 
found. The resonant angles are found to librate only when 
they are formed through convergent migration as expected. 
Note too that results obtained for two planet systems may 
also be relevant to systems of higher multiplicity if these can 
be regarded as being built up sequentially. 

Systems exhibiting second order commensurabilities are 
few in number and little studied in the context of migration 
models. However, further consideration may indicate that 
they could reveal information about their origin, so we have 
considered them in some detail. From the confirmed exo¬ 
planets, we list systems with period ratios that could poten¬ 
tially place them in second-order resonances in Table [T] All 
of these have a fractional deviation from a commensurabil¬ 
ity, < 6 x 10 -3 , which might allow them to be within it 
( see Section EH below). However, note that the fractional 
deviation for Kepler 365 is less than 6 x 10 -4 and for Kepler 
29 and Kepler 87 it is less than 4 x 10 -5 . 

The discovery of two planets eithe r in or close to 7:5 res - 
onance in HD 41248 was announced bv ljenkins et al.l d2013ll . 
This subsequently became controv ersial resulting in some 
uncertainty as to their status (see I Jenkins fe Tuomil I2014I . 
and references therein). Note that the innermost planet in 
Kepler 87 is a giant. Although the focus of this paper is on 
the Earth to Super-Earth mass range, this was included for 
completeness. 

The plan of the paper is as follows. We begin by dis¬ 
cussing the orbital properties of planets of the mass range 
of interest in a second order resonance. In Sections I21-I2H 
we provide a semi-analytic discussion, giving expressions for 
the eccentricities as a function of the ratio of circularization 
time to the migration time for the case when the planets 
undergo self-similar migration while maintaining the reso¬ 
nance. These are found to give results in good agreement 
with numerical simulations. In addition we show that these 
resonant configurations can occur when the eccentricities are 
small provided migration rates are sufficiently small. We go 
on to give a brief discussion of higher order resonances in 
Section EH 

In Sections EMU we give the details of the general 
setup of the numerical surveys. We then go on to present 
the outcomes of the surveys, as outlined above, in Sections 
141- 14.61 Finally, in Section El we summarise and discuss our 
results. 
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2 A SEMI-ANALYTIC SOLUTION FOR 
PLANETS IN SECOND-ORDER 
RESONANCES UNDERGOING MIGRATION 
AND ORBITAL CIRCULARIZATION 

Here we use the formalism of iMurrav fe Dermottl d 1999h . 
These authors give the equations for the orbital elements of 
the outer planet adopting a polar coordinate system with 
origin located at the primary mass as 

da,2 2 dR 
dt ri202 d\2 ’ 

de 2 \A -e\QR y/l - e|(l - yj\ - e|) OR 
dt ri2a%e2 dw2 n 2 a|e 2 8X2 ’ 

dvo 2 _ \A — e 2 <9R 
dt ri2a%e2 de2 ’ 

dA 2 2 <9R , y/l-eZ(l-y/l-e$) OR 

dt 71202 aa2 ri2a^e2 ae2 

where 02 is the semi-major axis, e 2 the eccentricity, VJ 2 the 
longitude of the periapsis and A 2 the mean longitude of the 
outer planet. Here and in the following sections, we iden¬ 
tify the two planets by applying a subscript 1 to denote the 
inner planet and a subscript 2 to denote the outer planet. 
For the outer planet R = RnGmi/a 2 , with mi being the 
mass of the inner planet and G the gravitational constant. 
Here Rd is the direct part of the disturbing function. The 
indirect part does not contribute significant terms in the 
discussion presented below. The equations governing the or¬ 
bital elements of the inner planet take the same form as o 
but with a 2 ,e 2 ,tU 2 , and A 2 replaced by ai,ei,z&i, and Ai 
respectively. In this case R is replaced by R_dGto 2 /o 2 . 

The quantity Rd may be developed as a Fourier series 
in A i and rut, i = 1,2. In order to discuss the second order 
resonance for which (p + 2 )ri 2 ~ pn\ where p is an integer 
^ 3, we retain terms containing Ai and A 2 only in the com¬ 
bination 8j = (p + 2)A2 — pAi. Thus the retained part of Rd 
takes the form 

Rd = FT cos (ndj — izu 1 + (* — 2n)vj2), (2) 

n,i 

where n is a positive integer or zero, i is a positive or negative 
integer or zero and the amplitudes FT are functions of ei, e 2 
and a = < 11 / 02 . For second order resonances, the lowest or¬ 
der amplitudes are homogeneous second degree polymomials 
in the eccentricities with coefficients that are functions of a. 
The forms of the FT c an be found correct t o fourt h order in 
the eccentricities from IMurrav fe Dermotl] (Il999l 'l. See Ap¬ 
pendix [A] below for more explanation. 

Note that the terms with n = 0 in © correspond to 
the standard secular terms which have to be retained in the 
discussion of second order resonances considered here. For 
second order resonances, the general angle 
ndj — izJi + (i — 2 n)w 2 occurring in 0 can be expressed as 
a linear combination of the two resonant angles 

0i = (p + 2)A2 — pAi — 2w\ and (3) 

03 = (p + 2)A2 — pAi —ZU 1 —ZJ 2 (4) 

For example 0 2 = 203 — 0i = (p + 2)A2 — pAi — 2vj2 and 
ZU 1 — ZU2 = 03 — 01. 

Following on from this it can be seen that a complete 
set of six equations for the state variables ai,et,i = 1,2, 


and 0 i, and 03 may be obtained by taking the first two 
equations of 0 and their counterparts for planet 1 , together 
with equations for 01 and 03 that may be obtained from the 
third and fourth equations of 0 and their counterparts for 
planet 1 . 


2.1 Migration and orbital circularization 


When a gaseous disc is present, the interaction of a planet 
with it will result in changes to its orbital energy and eccen¬ 
tricity. Thus if the orbital energy dissipation rate associated 
with the interaction of the disc with planet i is Di, the cor¬ 
responding rate of change of the semi-major axis is given 
by 


dai 2 0 ,■ 

~dF = ~ GM,rrii i 


(5) 


In addition, tidally induced orbital circularization causes e; 
to change at a rate given by 


da 

dt 


/ Tc,i , 


( 6 ) 


where r c ,i is the circularization time of planet i. To take 
account of interaction with the disc, the rates of change given 
by @ and © are added to the expressions for dai/dt and 
dei/dt given by the first two equations of 0 for planet 2 
and its counterpart for planet i. 


2.2 Self-similar migration 

We consider evolution of the system corresponding to res¬ 
onant self-similar migration in which 02 / 01 , ei, and e 2 are 
constant, the angle 03 is close to 0 and 0 i is close to n. 
Thus the semi-major axes decrease on account of migration 
torques while maintaining the same ratio. Solutions of this 
type have been considered by Papaloizou (2003) and Pa- 
paloizou & Szuszkiewicz (2005) for the case of first order 
resonances and as the generalisation to consider higher or¬ 
der resonances is straightforward in principle, the reader is 
referred to those papers for details omitted here. 


2.3 Conservation of angular momentum and 
energy 


The total angular momentum of the planetary system J = 
Ji + J 2 is given by 

Ji + J 2 = mi-\/GM*ai(l — e\) + m 2 \/GM*a 2 (l — e|) (7) 


and the total energy E is given by 


GM*mi 

2ai 


GM*m 2 

2<l2 


— m,2R 


(8) 


In the absence of migration and circularization these are 
both conserved. 

Following Nelson & Papaloizou (2002), for self-similar 
migration induced by interaction with the disc we have 


dJ _ j 1 dai 
dt 2ai dt 

^mi-^/GlVLai)! 



rri 2 yja 2 (\ - e|) 
7711-^/01(1 - el ) 

- el) ^ m2-^/GM*a2(l — 




(9) 




'T’mig, 2 
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where the total torque T = T\ + T 2 is the sum of the con¬ 
tributions of the migration torques acting on mi and m 2 . 
Here r m i g ,i is defined to be the migration rate of planet 
i. For an isolated planet with zero eccentricity we have 
d,Tli /dt = <$Tli /Tmig,i • 

Neglecting the small contribution of the interaction en¬ 
ergy, —m 2 R, the conservation of total energy gives 


dE_ 

dt 


GM,mi da\ 
2 af dt 


/l + m20i\ 

V mia 2 J 


= —D\ — D 2 = —D, 


( 10 ) 


where the total tidally induced orbital energy loss rate, D — 
D 1 + D 2 , is related to the torques acting on the planets and 
their circularization times through 

n _ GM*miei + GM*m 2 el n 2 T 2 mTi 
°i(!-ef)r c ,i o 2 (l - e|)r c , 2 yjl - e\ ' 

By eliminating dai/dt from © and (HOI) we can obtain a 
relationship between ei, e 2 , r c , Tmi g ,i and r m i g , 2 which with 
the help of m becomes 


1 f m,ia 2 e\ + m 2 aie| 
mim 2 ^^, 1(1 — el) r c , 2 (1 —e|) 

«2 7 V(! ~ 4) ~ a3 / 2 Vi 1 ~ e i) f _J__^ 

m 2 y / a 2 (l — e^) + mi ^/ai(l — ef) \ r mi g ,2 Tmi g ,i J 



We comment that this is general in that it depends only on 
the conservation laws and applies for any magnitude of ec¬ 
centricity or type of resonance. Note too that the right hand 
side of equation (fl2ll is proportional to l/r m i s , 2 — l/r m i 9 ,i 
which is the difference between the migration rates of the 
outer and inner planets. This should be positive, correspond¬ 
ing to convergent migration, if this equation is to be satisfied. 
However, self-similar migration with equilibrium eccentrici¬ 
ties was assumed, it may not apply when the eccentricities 
grow continuously. 


2.4 An expression for the eccentricity ratio 

We obtain an equation for the ratio of the eccentricities of 
the two planets by using © to find an expression for the 
rate of change of /3 = $3 — (f>i which takes the form 


dfi _ 

dt ^ 

n,i 

where the angle 1 /)" = n 6 j — ivj 1 + ( i — 2 n ) zu 2 
= ( 2 n - i) <j > 3 - (n - i) fa. 

For the self-similar migration we consider, <f> 1 is assumed to 
be very close to 7 r and <j >3 very close to 0. As the angles 
appearing in m are linear combinations of these angles, 
the right hand side is readily evaluated in this limit. The 
deviations of (pi and <p 3 from 7 r and 0 are small, vanishing 
in the limit of large migration time and changing no faster 
than on that time scale. Accordingly in that limit, d/3/dt 
may be neglected in m, which, as a is determined by the 
resonance condition, reduces to an algebraic equation for the 
eccentricity ratio ei/e 2 . For details about the determination 
of the F.J 1 for the second order resonances we consider below, 
see Appendix El 


i/l - e\m 2 OF? e l m i dF.f 


mafa 2 ei dei 


n 2 a2e 2 


de 2 


cos ip? (13) 


2.5 Resonant angles 

Use of (1121) and (1131) enables the determination of ei and 
e 2 in terms of the migration and circularization time scales. 
When these are comparable for both planets, the eccentric¬ 
ities will scale as the square root of the ratio of the char¬ 
acteristic circularisation time scale to the characteristic mi¬ 
gration time scale (e.g. Papaloizou 2003). The deviations of 
cf) 1 and (p 3 from n and 0 may then be determined from the 
equations for dei/dt obtained from ©. In the limit of small 
eccentricities, for second order resonances these are found to 
be on the order of M */(miniT c ,i) which can be small for suf¬ 
ficiently large evolution times. The corresponding quantity 
for first order resonances is M*/(mini % /T m i 9 ,;T c ,i) which is 
smaller by a factor ~ 1/ei. Accordingly, solutions with small 
eccentricities are obtainable at more rapid migration rates 
in that case. 


2.6 Behaviour at small eccentricities 

In the case of first order resonances, resonant angles may 
continue to librate even when the period ratios differ signif¬ 
icantly from strict commensurability, being driven away for 
example by tidal effects causing orbital circularization ( Pa¬ 
paloizou & Terquem 2010, Papaloizou 2011), even when the 
magnitude of the eccentricities approaches zero. However, 
because Rd is of higher order in the eccentricities, this does 
not occur for higher order resonances and the period ratio 
has to remain close to strict commensurability for the reso¬ 
nant angles to librate. For second order resonances, one can 
estimate the possible deviation, in the limit of zero eccentric¬ 
ities, by considering the terms proportional to either mi or 
m 2 in Equation CCD- When the resonant angles are constant, 
these terms are both equal to ( p + 2)n 2 — pni. For the 5:3, 
7:5 and 9:7 resonances appearing in simulations below, for 
mi = m 2 , we find ni/n 2 — (p+ 2 )/p ~ A(mi/M®)(M 0 /M*), 
where K = 3 x 10 _o , 3.5 x 10 —5 , and 4.2 x 10 -5 for p = 3, 5 
and 7 respectively. This indicates that ni/n 2 — (p + 2)/p 
1.5 x 10 -4 in the Super-Earth regime. 

However, departures from strict commensurability can 
increase at larger eccentricities as the magnitude of the res¬ 
onant part of the disturbing function increas es. Then from 
Equation (8.58) of iMurrav fc Dermottl ()l999l l , we estimate 
ni/n 2 — (p+2)/p ~ A e (ei/O.l)y / (mi/M®)(M 0 /M*), where 
K e = 1.55x 10 — 3 ,2.1 x 10” 3 , and 2.6 x 10 -3 for p = 3, 5, and 
7 respectively. This indicates a libration width more than an 
order of magnitude larger for eccentricities ~ 0 . 1 . 


2.7 Higher order resonances 

Although the main focus of this section has been on sec¬ 
ond order resonances, we have observed the occurrence of 
even higher order resonances in a few of our simulations. 
We comment that the discussion of second order resonances 
can be generalised to apply to higher order resonances. For 
a fc-th order resonance the angles (pi and (j >3 are modified to 
become 

4>i = (p + k )\ 2 — pAi — kzu 1 and (14) 

<t>3 = (P + k )\ 2 - p\l - (fc — l)tUl — VJ 2 . (15) 

The formalism developed from © and © can be applied as 
before with all angles now expressed as linear combinations 
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of the new angles. However, the form of the FJ 1 , for n ^ 0 is 
at least of order k, whereas when n = 0 corresponding to the 
secular terms, as in the case of second order resonances there 
are second order terms. The above means that when the ec¬ 
centricity is not too small, the resonant angle deviations are 
expected to be ~ [M* / {rrnriiTmig , i)](Tmi g ,i/T c ,i ) k/2 , whereas 
for small eccentricities the system is expected to be domi¬ 
nated by secular interaction. 


3 NUMERICAL SIMULATIONS 


We have performed numerical simulations of a pair of copla- 
nar migrating low-mass pla nets using the N-body code RE¬ 
BOUND dRein fe Liul 1201^1 . Additional acceleration terms 
that model the effects of orbital migration and eccentricity 
damping that can be assumed to result from the interaction 
of the planets with a gaseous protoplan etary disc (see eg. 
iLee fc Pealdl2002 ; ISnellgrove et alJboQl .l are included. The 
resulting acceleration r of a planet is given by 


2rr rxj 2r r r ■ <p 

"F I 12 S 

Tc Tmig |l*| Tc Tmig 


(16) 


where Tmig corresponds to the timescale for migration 
adopted above and similarly t c to the timescale for orbital 
circularization. Note that we have dropped the subscript, i, 
that indicates the planet from these particular quantities. 
The force per mass arising from gravitational interactions 
with other bodies is f g , j is the angular momentum per unit 
mass, and r and f> are the unit vectors in the radial and 
azimuthal directions for polar coordinates with origin at the 
central star. 

Depending on the detailed properties of the disc, plan¬ 
ets in the mass range we consider can either undergo Type I 
migration in a linear regime, or be in a nonlinear regime 
for which part ial gap formation and/or wak e-planet inter¬ 
actions occur (I Baruteau fe Papaloizoull2013l l. Furthermore 
in the inner regions close to the star where truncation due 
to a magnetic field may occur, the disc structure may be 
strongly affected and strong outward propagating density 
waves may be present. Planetary migration and circulariza¬ 
tion rates depend on the radial form of the state variables 
and the equation of state in the protoplanetary disc (e.g. 
I Paardekooper et al.l 120101 ). These time dependent quanti¬ 
ties determine both the magnitude and sign of the migration 
rate which thus have considerable uncertainty over the disc 
lifetime which is characteristically 10 6 ~ 7 yr. 

In order to perform simulations on this characteristic 
time scale, we have considered a very simple disc model for 
which Tmig,i and T c ,i, i = 1,2, are constants exterior to some 
cavity radius with 1/ Tmig, i and 1 /r c ,i vanishing inside it in 
some cases. In other cases circularization was allowed to con¬ 
tinue for one or both planets while inside the cavity in order 
to model the situation where that continues after migration 
ceases. This can occur when the mean disc surfa ce density 
increases outwards fe.g. l Paardekooper et all201Ch . We have 
then surveyed a large range of possible values for T m i g ,i and 

Tc,i/Tmig,i - 


Survey label 

m\ 

mF 

m 2 

a i 

1AU 

&2 

1AU 

Extended 

N S ec 

lAU-l-4.21 

1 

4 

0.6 

1 


15 

1AU 1 4 32 

1 

4 

0.7 

1 

y 

10 

1AU 1 2 21 

1 

2 

0.6 

1 

y 

10 

1AU 1 2 32 

1 

2 

0.7 

1 


2 

1AU 2 1 21 

1 

2 

0.6 

1 


5 

1AU 2 1 32 

1 

2 

0.7 

1 


3 

5AU.1.4.21 

1 

4 

3.0 

5 


6 

5AU.1.4.32 

1 

4 

3.5 

5 


4 

5AU 1 2 21 

1 

2 

3.0 

5 


3 

5AU 1 2 32 

1 

2 

3.5 

5 


6 

5AU 2 1 21 

1 

2 

3.0 

5 


4 

5AU_2_1_32 

1 

2 

3.5 

5 
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Table 2. Parameters for the surveys. The first column gives the 
survey label. Successive columns then give from left to right, the 
masses of the inner and outer planets in earth masses, the ini¬ 
tial semi-major axes of the inner and outer planet in AU, and y 
indicating that the survey was extended with one or both plan¬ 
ets undergoing circularization inside the cavity. The final column 
gives the number of second order resonances present at the end 
of all the runs constituting the survey before any extension. 


3.1 Numerical surveys 

Each of the surveys that we present in this paper includes 
20 x 20 runs each of which is for a particular choice of T m i g ,i 
and Tc,i / Tmi g ,i ■ These quantities are identified by the pa¬ 
rameters Ni and N% which are generate values for these 
timescales for a planet of mass m,i through the equivalent 
quantities 

Tmig = 1.62^ x 10 3 x —^ = r m o—^-yr and 
rrn rm 

TgjTmig = 1.3^ X 3 X 10" 4 . (17) 

The values for Ni and N 2 adopted were zero or integers 
in the interval [0,19]. This leads to a migration timescale 
in the range (10 3 — 9.6 x 10 6 )(M®/mi) yr and the ratio 
Tc/Tmig in the range (3 x 10 -4 — 4.4 x 10 -2 ). The ratio 
T c /Tmig is expected to be ~ 2.5 x 10~ 3 for a disc with aspect 
rati o ^0.05 from a discussion of d isc-planet interactions (see 
e.g. IPapaloizou fe Larwoodl l2000l l. This lies in the middle 
of the considered range for one earth mass. We consider 
surveys for which the outer planet starts on a circular orbit 
at either 1 AU or 5 AU with the inner planet close to but 
outside either the 2:1 or the 3:2 resonance. In addition, we 
incorporate an inner cavity of radius 0.08 AU inside which 
migration torques do not operate. But note that our results 
can be scaled to other values using the scaling of the units 
of space and time available for systems governed by gravity 
(see Section ^. 2 II . Parameters associated with the surveys we 
carried out are listed in Table [5] 

For the first set of surveys, when a planet enters the 
inner cavity, the migration and circularization forces are 
switched off and it is only the gravitational forces between it 
and the other planet and the central star which determines 
the subsequent evolution. Runs are found to undergo little 
change in semi-major axis or period ratio once both planets 
are inside the cavity. As these are the main quantities of 
interest, the evolution can be regarded as being complete at 
this stage. Accordingly we stopped a run once both planets 
were interior to 0.075 AU. A simulation time of 10 7 yr was 
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found to be adequate to ensure the completion of all our 
runs. 

For the second set of surveys, a survey of the type de¬ 
scribed above is extended from the point at which both plan¬ 
ets are in the cavity with circularization forces continuing to 
apply only for the outer planet but with migration deacti¬ 
vated for both planets. The device of enabling circularization 
to continue for the outer p lanet once convergent mi g ration 
has ceased was adopted by I Baruteau fc Papaloizoul (l2013l l 
to simulate the effects of wake-planet interactions. It can be 
regarded as equivalent to a situation where there is a zero 
torque condition in the vicinity of the cavity boundary with 
circularization continuing to operate there. Such a cavity 
may be assumed to exist in practice as a consequence of a 
central magnetosphere or some other structural feature of 
the disc. 

3.2 Initial configuration 

We considered planets with masses of either 1, 2 or 4 earth 
masses (M®). The planets were initialised on circular copla- 
nar orbits with the outer one being at 1 AU or 5 AU. The 
inner planet is initialised on an orbit which is exterior to 
either 2:1 or 3:2 mean motion resonance. The initial period 
ratio is P 1 /P 2 = 2.15 in the former case and P 1 /P 2 = 1.7 in 
the latter case. Note that the system is invariant to multi¬ 
plying all lengths by a factor f s and times, including Tmi gi i 
and r c ,i , i = 1,2, by a factor (/ s ) 3 ^ 2 - These results can be 
extended to apply to cases with a different cavity radius. 

4 RESULTS 

4.1 General features 

In all the unextended surveys listed in Tabled the outcome 
is that either the planets are in resonance or they collide. 
First and higher order resonances are possible but the lat¬ 
ter are infrequent. We note that although planets attaining 
first order resonances may u ndergo an instability that cause s 
the resonance to be broken dGoldreich fc Schlichtinall2014ll . 
another resonance with smaller period ratio is quickly at¬ 
tained. Thus most of the time is spent in some resonance 
unless the migration is fast enough to cause the planets to 
collide. 

All our runs are stopped at that point. To illustrate 
a typical case, we describe the evolution of the period 
ratio and eccentricity for a run taken from the survey, 
1AU_1_4_21, for which the complete set of results are il¬ 
lustrated in Fig. [I] The final period ratio is indicated in 
the eighth row and eighth column in the upper leftmost 
panel. The results are illustrated in Fig. [2] We recall that 
the outer planet mass was 4 M® and the inner planert mass 
1 M®. Our migration prescription ensures convergent mi- 
gration that first leads to trappin g in the 3:2 resonance (e.g. 
IPapaloizou fe Szuszkiewiczl l2005l l. However, as the migra¬ 
tion continues this is broken and a 4:3 resonance attained. 
That in turn is broken leading to a 5:4 resonance which re¬ 
mains until the pleanets both enter the cavity. Note that 
towards and until just before the end of the run, the ec¬ 
centricities increase. This is a common feature of the runs 
presented here, being due to the loss of the damping effect 


of the circularization of the inner planet while it is tem¬ 
porarily driven inwards by the outer planet. Note that the 
transition from reson ance to resonance seen here was a lso 
seen in simulations bv IPapaloizou fe Szuszkiewicd (l2009l l. 


4.2 Results of the surveys with 

(mi,m 2 ) = (1,4) M® 

In Fig. □ the final period ratios for planet pairs initialised 
on orbits close to the 2:1 resonance are shown in the left set 
of four panels and close to the 3:2 resonance in the right set 
of four panels. For each set of four the upper and lower left 
panels give the result when the outer planet starts at 1 AU 
and the upper and lower right panels are for the outer planet 
starting at 5 AU. Apart from a few exceptions associated 
with the planets scattering as they enter the cavity, runs 
end with either the planets close to a resonance, or with a 
collision. Most of the resonances are first order and readily 
identified in Fig. |T]and similar figures below with the help 
of the colour bars (e.g. for surveys starting close to the 2:1 
commensurability illustrated in Fig.Q] a filled dark red circle 
in the lower leftmost panels corresponds to 2:1 resonance, a 
filled purple circle in those two panels to the 3:2 resonance, 
a filled orange circle in the upper leftmost two panels to 4:3 
resonance and a filled light green circle in the latter panels 
to 5:4 resonance). 

In general planetary collisions are expected to occur for 
very fast migration (small migration time). However, the 
rate of migration required depends on the rate of orbital 
circularization, decreasing as the rate of circularization de¬ 
creases. This is because weak circularization allows for the 
development of highly eccentric planetary orbits, making a 
collision through orbit crossing more probable. This trend 
can be seen in Figure [T] (and other surveys presented be¬ 
low). As migration times increase, collisions tend to require 
larger circularization times in order to occur. For exam¬ 
ple, from the left group of four panels in Fig. [X] we can 
quantify the migration and circularization time limits for 
collision. In the left panels with the outer planet start¬ 
ing at 1 AU, collisions can occur for r c /r m i g ^ 0.01 pro¬ 
vided that TmO ^ 10 4 yr. This corresponds to the quan¬ 
tity, l/{l/rmig,2 - l/Tmi g ,i) = T m o/3, which estimates the 
convergence time through migration of the orbits of the 
outer and inner planets (see equations m and 113), being 
^3x 10 3 yr. However, when r c /r m i g ~ 0.1, collision can oc¬ 
cur for Tmo ~ 10 6 yr. Note that these times approximately 
increase by a factor 5 3 ^ 2 for the right group of panels with 
the outer planet starting at 5 AU. 

To decide whether a run ended with a relatively rare 
second order commensurability, we adopted the same crite¬ 
rion as above that fractional deviation from commensurabil¬ 
ity should be < 6 x 1CU 3 . This was applied to all surveys. 
These were indicated in Fig.|T]and other similar figures by a 
cross. Thus for surveys with the outer planet starting close 
to the 2:1 resonance, a 5:3 resonance was indicated by an 
olive cross in the lower panel that is second from the left, 
a 7:5 resonance by a burgundy cross in the two upper left¬ 
most panels and a 9:7 resonance by a yellow cross in the 
upper panel that is second from the left. These are discussed 
in more detail in Sections and I4.4.TI below. The survey 
1AU_1_4_21 had also had runs that were associated with 
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Figure 1. Final period ratios for planet pairs initialised on orbits close to 2:1 (left group of four panels) and 3:2 (right group of four 
panels) with the inner planet mass being 1 Mq and the outer mass being 4 Mj. For each group of four panels, the upper and lower 
left panels give results for the outer planet starting at 1 AU and the upper and lower right panels give results for the outer planet 
starting at 5 AU. These respectively correspond to the surveys 1AU_1_4_21, 5AU_1_4_21, 1AU_1_4_32, and 5AU_1_4_32. Period ratios 
are indicated in associated colour bars such that the larger values up to 2 are indicated in the lower panels and the smaller values down to 
1 are indicated in the upper panels. Second order resonances are indicated by crosses (see text) and the purple squares indicate collisions 
between the planets. Note that in this and other similar figures, the migration time plotted is r m o and the ratio of circularization time 
to migration time plotted is T c /T m i g which are both independent of planet mass (see Equations 117H . For the simulations shown in the 
above panels the quantity, l/(l/r m i 9i 2 — 1 ), which estimates the initial convergence time for the outer and inner orbits is equal 

to Tmo/3. 




Time [yr] Time [yr] 

Figure 2. The evolution of the period ratio (left panel) and eccentricity (right panel) for a run contributing to the survey 1AU_1_4_21 
(see text). The vertical line in this and similar figures indicates the point at which both planets have entered the cavity. 


higher order resonances with period ratios close to 10:7, 8:5 
and 11:6. These are discussed in more detail in Section ISO 


4.3 Results of the surveys with 

(mi, m 2 ) = (1, 2) M® and (mi, m 2 ) = (2,1) M® 

These surveys are carried out and the results are presented in 
the same manner as the surveys with (mi, m 2 ) = (1,4) M®, 
the only change being that they are carried out with different 
planet masses. The qualitative description of the results is 
found to be the same. 
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Figure 3. Final period ratios for planet pairs initialised on orbits close to 2:1 (left group of four panels) and 3:2 (right group of four 
panels) with the inner planet mass being 1 Mq and the outer mass being 2 Mg. For each group of four panels, the upper and lower 
left panels give results for the outer planet starting at 1 AU and the upper and lower right panels give results for the outer planet 
starting at 5 AU. These respectively correspond to the surveys 1AU_1_2_21, 5AU_1_2_21, 1AU_1_2_32, and 5AU_1_2_32. Period ratios 
are indicated in the associated colour bars such that the larger values up to 2 are indicated in the lower panels and the smaller values 
down to 1 are indicated in the upper panels. Second order resonances are indicated by crosses (see text) and the purple squares indicate 
collisions between the planets. For these cases, the quantity l/(l/r m i gi 2 — l/m m i 9 ,i), which estimates the initial convergence time for 
the outer and inner orbits is r m o, which is the migration time plotted. 



time [yr] time [yr] 


Figure 4. The evolution of the period ratio (left panel) and eccentricity (right panel) for a run contributing to the survey 5AU_l_2_32(see 
text). In this case the planets undergo orbit crossing and scatter causing the inner and outer planets to be interchanged. 


The results for the surveys with (mi,m 2 ) = (1,2) M® 
are shown in Fig. [3] These surveys produce a similar number 
of collisions, range of first order resonances and second order 
resonances as the previous ones. However, they differ in that 
a significant number of 2:1 commensurabilites are formed 
at high migration rates and relatively low circularization 
rates. These are seen in the second and fourth of the lower 
panels counted from the left of Fig. [3] They are found to be 
associated with a scattering and interchange of the planets. 


Figure [3] shows the evolution for a run taken from the survey 
5AU_1_2_32. It is indicated in the ninth row from bottom 
and rightmost column of the rightmost lower panel of Fig. 
[21 In this run the planets undergo orbit crossing and scatter 
causing the inner and outer planets to be interchanged and 
a large period ratio ~ 14 attained. Subsequently convergent 
migration resumes and a 2:1 commensurability is formed. 

Results for the surveys with (mi, m 2 ) = (2,1) M® are 
shown in Fig. [5] In this case the general form of the or- 
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Figure 5. Final period ratios for planet pairs initialised on orbits close to 2:1 (left group of four panels) and 3:2 (right group of four 
panels) with the inner planet mass being 2 Mq and the outer mass being 1 Mj. For each group of four panels, the upper and lower 
left panels give results for the outer planet starting at 1 AU and the upper and lower right panels give results for the outer planet 
starting at 5 AU. These respectively correspond to the surveys 1AU_2_1_21, 5AU_2_1_21, 1AU_2_1_32, and 5AU_2_1_32. Period ratios 
are indicated in associated colour bars such that the larger values up to 2 are indicated in the lower panels and the smaller values down to 
1 are indicated in the upper panels. Second order resonances are indicated by crosses (see text) and the purple squares indicate collisions 
between the planets. The quantity l/(l/r m i ffi i — 1 /T mv , y 2 ), which for these cases estimates the initial divergence time for the outer and 
inner orbits is r m o, which is the migration time plotted. 


bital evolution differs from the previous cases. As the inner 
planet is more massive than the outer planet, the initial in¬ 
ward migration is divergent. Only when the inner planet 
has reached the inner cavity while the outer planet is still 
migrating can their semi-major axes start to approach each 
other and their period ratio starts to decrease so that a 
resonance can form. The nature of this evolution allows a 
significantly larger number of 2:1 comensurabilities to result 
from systems with larger migration rates, as compared to 
the previous surveys. 

Note that the simulations with the outer planet starting 
at 5 AU, would lead to the same results as in the case of 
the outer planet starting at 1 AU if the timescales were 
multiplied by a factor of 5 3 ^ 2 and the inner cavity radius 
was also increased by a factor of five. However, because the 
inner cavity radius is in fact the same in both cases, the 
simulations with the outer planet starting at 5 AU map into 
simulations starting at 1 AU with a reduced inner cavity 
radius. This results in a simulation that is allowed to evolve 
for longer. This has the consequence that the planets end 
up closer together or with smaller period ratios than would 
occur if the simple scaling was assumed. This is apparent for 
example when the left and right panels in the leftmost group 
of four panels in Fig. [5] are compared. However, the simple 
scaling is relatively more closely adhered to when collisions 
are considered because the inner cavity plays less of a role 


in those cases (see discussion of the results shown in Figs[l] 
in section POP . 

In Figurc[S]we show the orbital evolution of a run taken 
from the survey 1AU_2_1_32. It corresponds to the run indi¬ 
cated in the eleventh row (from bottom) and eleventh col¬ 
umn (from left) of the third lower panel from the left in Fig. 
[5] During the first part of the simulation, the inner planet 
migrates much faster inwards than the outer planet leading 
to a significant increase of their period ratio. This growth 
suddenly comes to a stop when the inner planet has reached 
the cavity and its migration switches off. From this time on, 
the period ratio decreases until a persisting 3:2 mean-motion 
resonance forms. 

In addition to producing more 2:1 commensurabilities 
significantly fewer collisions were found in these surveys as 
compared to the previously discussed ones. However, a com¬ 
parable number of second order resonances were found. 

4.4 Occurence of second-order resonances 

In Table [2] we show the number of second-order resonances 
found in all of the surveys discussed above. We detected 5:3, 
7:5 and 9:7 resonances. These are indicated in the figures 
illustrating the final period ratios. 

From Table [2] there were 72 occurrences out of a to¬ 
tal of 4800 simulations giving an incidence rate of 1.5%. 
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Figure 6. Evolution of the semi-major axis and the period ratio for a run taken from the survey 1AU_2_1_32 for which the inner planet 
is more massive (see text). 


Resonance 

ei 

e2/ei 

r 

5:3 

0.13 

0.07 

0.65 

0.571 

0.571 

0.557 

7:5 

0.087 

0.575 

0.504 

9:7 

0.075 

0.03 

0.560 

0.5 

0.456 

0.481 


Table 3. Comparison of numerical (third column ) and semi- 
analytically determined (fourth column) eccentricity ratios. The 
resonance and value of the eccentricity of the inner planet are 
indicated in the first and second column, respectively. 

Although, these occurrences are rare, as they may contain 
useful information concerning their origin, we now examine 
some of them in more detail. 



0.0e+00 5.0e+05 1.0e+06 1.5e+06 

Time [yr] 

Figure 8. The period ratio evolution for the run taken from 
5AU_1_4_32 that shows a scattering followed by the attainment 
of a 5:3 resonance (see text). 


4-4-1 Detailed study of second-order resonances 

We now examine some runs which ended with the planets 
in a second order resonance in more detail. For these the 
inner planet mass was 1 M® and the outer planet mass was 
2 M®. The outer planet started at 1 AU just outside the 
2:1 resonance with the inner planet. They form part of the 
survey 1AU_1_2_21 for which results are illustrated in the 
panels furthest to the left of Fig. [3] We consider the unique 
runs in this survey that ended with the planets in 7:5 and 
9:7 resonances and the run corresponding to the second to 
last entry in the fifth row from the top of the lower panel 
furthest to the left. The latter (case A) is one of several that 
ended in 5:3 resonance. It had T m ig = 1.5 x 10 6 (M®/mi) yr 
and r c /Tmig = 0.03. For the 7:5 resonance case (case B), 
Tmig = 1.2 x 10 s (M®/m;)) yr and TcfTmig = 0.015 and 
for the 9:7 case (case C) r m i 9 = 5 x 10 4 (M®/mi)) yr and 

T c /Tmig = 0.012. 

The evolution of the semi-major axes, period ratios, ec¬ 
centricities and resonant angles for these cases is illustrated 
in Fig. □ At early times a first order resonance is attained, 
2:1 in case A, 3:2 in case B and 4:3 in case C. This be¬ 
comes unstable and the second order resonance is entered 
and maintained until the end of the run. 0 i, 02 and 03 are 
the resonant angles as defined in Eqs. © - ©. For all three 
simulations, when the second-order resonance occurs, (f> 1 and 
02 are around n while 0 3 is librating around 0 and 27r. We 
recall that 0 J 2 — vji is the angle between the apsidal lines 
which is found to librate around n. 


Averaging out short period oscillations, constant equi¬ 
librium eccentricities, as indicated in Section are at¬ 
tained, though this phase is short lived in case B. The ec- 
centricties start to jump just after the inner planet starts 
to enter the cavity. This is because when it no longer inter¬ 
sects the disc, migration and circularization then only oper¬ 
ate on the outer planet, At later times, final constant values 
of the eccentricities are attained that are larger than those 
applying when both planets were affected by migration and 
circularization. 

We have tested the applicability of Equation (HI) for 
the eccentricity ratio and Equation HI which, with the help 
of Appendix m enables the calculation of the magnitues 
of the eccentricities during phases when the system under¬ 
goes self-similar migration, which we here identify with the 
phases when the eccentricities are approximately constant 
with superposed small amplitude short period oscillations. 
Equation m then enables the calculation of the eccentrici¬ 
ties as a function of the migration and circularization times. 
In addition comparison of the numerically determined ec¬ 
centricity ratios with those obtained from the semi-analytic 
procedure is given in Table [3j It is seen that there is good 
agreement, especially for the smaller eccentricties. In these 
cases adequate calculations can be performed taking only 
second order terms into account (see Appendix [©. The ab¬ 
solute values of the eccentrities are also found to agree with 
those implied by m- 

In Fig. HI there are two green crosses corresponding to 
5:3 resonances that can be seen in in the rightmost lower 
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Figure 7. Semi-major axis (first column), period ratio (second column), eccentricity (third column) and resonant angles <j> 1 , <j >2 and 03 
as well as m 2 — m\ (fourth column) for runs that attained the commensurabilities 5:3 (upper panels), 7:5 (middle) and 9:7 (lower panels). 
The point at which both planets have entered the inner cavity is indicated by the vertical line. 


panel showing the survey with the outer planet starting at 5 
AU. In these cases the relatively rapid convergent migration 
rates result in orbit crossing leading to the planets changing 
places. The original inner planet is scattered outwards such 
that the system attains large period ratios > 14. Convergent 
migration resumes after the now more massive inner planet 
enters the cavity resulting in the attainment of a 5:3 reso¬ 
nance. The period ratio as a function of time for the case 
with the shortest migration time is shown in Fig. [8] 


4.5 Runs indicating higher order resonances 

To illustrate the potential presence of higher order reso¬ 
nances we illustrate the evolution of runs taken from the 
survey 1AU_1_4_21. These were cases with final period ra¬ 
tios 1.83 (case D), 1.6 (case E) and 1.427 (case F). Case D 
corresponds to the run shown in the fourth row (from top) 
and third column from the end in the lower leftmost panel 
of Fig. [T] Case E corresponds to the run shown in the top 
row and fourth column from the end in the lower leftmost 
panel of Fig.Q] Case F corresponds to the run shown in the 
fourth row and second column from the end in the upper 
leftmost panel of Fig. [T] The time independent evolution of 
the semi-major axes, eccentricities and angle between the 
apsidal lines for these cases is illustrated in Fig. [9] Each of 
these r uns initially attains a 2: 1 reson ance which is not sus¬ 
tained (iGoldreich fc Schlichtinell2014l ). The systems leave it 
and attain period ratios close to 11:6, 8:5 and 10:7 in cases 
D, E and F respectively. The angle between the apsidal lines 
oscillates between 0 and 7r before settling on 7r in cases D 
and E and 0 in case F after the planets have entered the 
cavity and forces arising from migration and circularization 
were removed. Note that the eccentricities are large in these 
runs attaining ~ 0.45 in case D. 


To investigate the subsequent evolution of the period 
ratios that could occur if additional dissipative forces were 
imposed we performed an additional set of runs for which 
circularization was allowed to continue to operate on the 
outer planet after the cavity was entered, The evolution of 
these runs is illustrated in Fig. 1101 It will be seen that the 
eccentricities decrease while changes to the period ratio are 
on average very small. The inner planet has the larger ec¬ 
centricity while dissipation occurs only for the outer planet. 
The eccentricity of the inner planet decreases through secu¬ 
lar transfer to the outer planet which is indicated by main¬ 
tenance of the angle between the apsidal lines at 0 or 7r. In 
case D the latter fails once both eccentricities become very 
small and the angle loses its significance. In all cases the 
evolution ceases as the eccentricity of the outer planet ap¬ 
proaches zero. During the later stages of evolution, energy is 
removed from the orbit of the outer planet, but not that of 
the inner planet, causing the period ratio to decrease very 
slightly. 


4.6 Continuation of runs with orbital 

circularization but with no migration torques 

In this section, we describe some surveys of the second type 
that were continued after the planets have reached the inner 
cavity. We begin by considering extensions of the surveys 
1AU_1_2_21 and 1AU_1_4_32 the results of which are shown 
in Figs [3] and [l] 

For the extensions to the runs, as described above we 
allow circularization to continue only for the outer planet 
while migration is deactivated for both planets interior to 
the cavity. Each run was then continued for a time which is 
the shorter of either 5 x tmax or 2 x 10 1 yr, with t max being 
the duration of the unextended run. 
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Figure 9. Semi-major axis (first column), period ratio (second column), eccentricity (third column) and the angle between the apsidal 
lines — pi (fourth column) for runs with ultimate period ratios near 1.83 (upper panels), 1.6 (middle panels) and 1.427 (lower panels) 
that were carried out in the survey 1AU_1_4_21. 
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Figure 10. As in Fig. [9] but in this case after point where both planets enter the cavity, the runs are extended with migration forces 
switched off while circularization forces continue to act only on the outer planet. 


Figure [lT] shows the final period ratio at start of this 
continuation on the left panel and the final period ratio after 
continuing the simulations. By making a comparison one can 
see that the period ratios are increased. This is expected 
because during this evolution, the planetary systems have 
energy removed at constant angular momentum, and the 
theoretical expectation is that they should separate with 


perio d ratio increasing as a function of time fe.g. l Papaloizou! 

ImlO). 

This is especially the case for small values of the ratio 
of the circularization time to the migration time (see the left 
hand sides of the panels). For example inspection of the two 
lower leftmost panels of Fig. [TT] shows systems that were in 
3:2 resonance increasing their period ratios in to the 1.52 - 
1.56 domain for both surveys illustrated. However, changes 
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Figure 11. Final period ratios after extended runs for planet pairs initialised on orbits close to 2:1 (left hand set of four panels) and 
3:2 (right hand set of four panels) are shown in the upper and lower right hand panels of each group of four. The upper and lower left 
hand panels for each set show the period ratios at the commencement of the extension of the runs ( the result of the unextended surveys 
1AU_1_2_21 and 1AU_1_4_32) for comparison. The inner planet mass is 1 M.j and the outer mass 2 Mj in the first case and 4 Mj in 
the second. The outer planet starts at 1 AU for all runs. The purple squares again show simulations with collisions between the planets. 


are much smaller for longer circularization times. In addi¬ 
tion the period ratio range 1.7 — 1.8 is empty. This can be 
rectified if shorter circularization times are used so induc¬ 
ing a more rapid separation of the planetary systems. We 
no te the p henomenon of wake-pla net interaction discussed 
by I Baruteau fe Papaloizoul (i2013l ) that could lead to the 
same effect and was modelled through the introduction of 
orbital circularization. 

To illustrate how the results depend on the amount 
of circularization (or wake-planet interaction) introduced 
we have perfomed a simplified illustrative calculation by 
extending the continuation of 1AU_1_4_32 described above 
such that each run is extended for a further 0.1 tmax, where 
tmax is the duration of the run in the unextended survey, 
during which the circularization ti me for each planet taken 
to be 20 orbits, as was adopted in i Baruteau fe Papaloizoul 
d2013l) . Results are shown in Fig. 1121 It will be seen that 
period ratios throughout the range 1.7 — 2 are obtained. 

The criteria for a run to have terminated in a second or¬ 
der resonance were applied in the same way as for the other 
surveys and more are found after this doubly extended sur¬ 
vey completes. However, as in this case they are formed from 
divergent migration, capture into libration is not expected. 
This is indeed found to be the case. Thus runs of the type il¬ 
lustrated in Fig. 0 which end with resonant angles librating, 
do not occur when the near second order commensurabili- 
ties are obtained as a result of divergent migration produced 
through the action of efficient orbital circularization. 


5 SUMMARY AND CONCLUSIONS 

In order to investigate the possible consequences for poten¬ 
tial architectures of low mass planetary systems, we have 
performed numerical simulations of a pair of coplanar mi¬ 
grating low-mass planets in the mass range (1 — 4) M®. In 
the surveys we conducted, the planets were initialised on 
circular coplanar orbits with the outer one being at 1 AU 
or 5 AU. Outside an inner cavity of radius 0.08 AU, the 
orbital evolution was considered for a constant migration 
time, Tmig, in the range (10 3 — 10 7 )(M®/mi) yr and con¬ 
stant ratio of circularization time to migration time r c /r m i 9 
in the range 3 x 10 -4 — 4.4 x 10“ 2 . These could be viewed 
as possibly arising from interaction with a protoplanetary 
disc. The inner planet was initialised on an orbit which was 
close to and exterior to either the 2:1 or the 3:2 mean mo¬ 
tion resonance. Note that the length and time scales we used 
can be scaled to other values using the standard scalings for 
systems governed by gravity (see Section [3.21) . 

In order to provide a simple model of phenomena pro¬ 
ducing an inner cavity boundary to a protoplanetary disc, 
such as rapid changes in disc viscosity or interaction with 
the stellar magnetic field, we considered migration to termi¬ 
nate inside the cavity but performed some extended surveys 
for which circularization was allowed to continue for one 
or both planets. We remark that the latter set up could be 
viewed as providing a heuristic modelling of phenomena such 
as wake-pl anet interactions, or planet-den sity wave interac¬ 
tions (e.g. I Baruteau fe Papaloizoul |2013j) or even residual 
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Figure 12. Final period ratios after the runs in the survey ex¬ 
tended from 1AU_1_4_32 that are shown in Fig. Ill I were further 
extended for an additional time, 0.1 tmax, with a circularization 
time of 20 orbits for both planets, are shown in the right hand 
panels. The upper and lower left hand panels show the period ra¬ 
tios at the commencement of the first extension of the runs (the 
result of the unextended survey 1AU_1_4_32). The purple squares 
show simulations with collisions between the planets. 


time. The ratio of the eccentricities was subsequently eval¬ 
uated correctly to fourth order in the eccentricities for sev¬ 
eral cases and was found to be in reasonable agreement with 
simulations. In a very small number of cases, period ratios 
corresponding to even higher order resonances were seen. 
These were discussed in Sections nm and HU 

In order to obtain an extended range of final period 
ratios, it was necessary to extend the simulations so that 
the planet (s) underwent a period of circularization without 
migration. The important consequence of this is that there 
is a period of effective divergent migration that causes the 
period ratios to diverge. Depending on its effectiveness, the 
whole range of period ratios in (1,2) could be obtained (see 
Section [4.611 . Thus, in order to determine final outcomes of 
planetary systems, the study of their interaction with the 
inner parts of a protoplanetary disc is crucial. A better un¬ 
derstanding of the detailed inner disc structure, the effects 
of any density waves launched by an interaction with the 
central star and the interaction of any residual planetesi- 
mal disc will be required. These issues will form the focus of 
future studies. 
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planetesimal migration (e.g. iMinton fe Levisonl 120141 ') that 
can halt or turn around convergent migration of the plan¬ 
ets. We remark that the work here is carried out in order to 
characterise the possible outcomes for migrating two planet 
systems under very simple assumptions about forces lead¬ 
ing to migration and circularization. It does not constitute 
population synthesis studies. 

Numerical simulations that terminated without an ex¬ 
tended period of circularization resulted for the most part 
in, either collisions between the planets or, the formation 
of systems close to resonance. These were mostly first order 
with second order resonances occurring in a few % of cases. 
Those were found to occur for modest eccentricities <~ 0.1 
and all resonant angles librated (see Section [OH]). This is in 
contrast to the situation where the near commensurability 
resulted from divergent migration of the type encountered 
in the extended surveys. In those cases resonant trapping 
does not occur and the resonant angles circulate. Thus the 
behaviour of the resonant angles in systems near a second 
order commensurability can indicate the form of any migra¬ 
tion that led to them. 

An analytic description of planets migrating in second 
order resonances was given in Sections [2] - 12.61 There we 
derived expressions for the ratio of the eccentricities of the 
planets and the absolute values of the eccentricities as a 
function of the ratio of circularization time to migration 
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APPENDIX A: THE DIRECT PART OF THE 
DISTURBING FUNCTION FOR SECOND 
ORDER RESONANCES 

In order to discuss the (p + 2) : p resonance, the quantity 
Rd is written in the compact general form 

Rd = Fcos (ndj — i-oji + (i — 2 n)vj 2 ), (Al) 

n,i 

where 6 j = (p + 2) A 2 — pAi. Here the index of summation n 
is a non negative integer and the index of summation, i, is a 
positive or negative integer or zero. In addition terms that 
turn out to have the same cosine factor are collected together 
in the expressions given below. We take into account terms 
up to fourth order in the eccentricity such that Rd takes 
the explicit form 

Rd = Fq + Fi° cos(u72 — tui) + F° cos 2(u72 — roi) + 

F 2 cos (9j — 2tui) + Fi cos ( 8 j — m 2 — roi) + 

Fq COS (9j — 2 tU 2 ) + F 3 1 COS (0j + tU 2 — 3zU\) + 

Flj cos ( 8 j + U 7 i — 3 zu 2 ) + F 4 cos(20j — 4zui) + 

F 3 cos(20j — 3zui — VJ 2 ) + F| cos(2 9j — 2 -cji — 2 zu 2 ) + 

Fi cos(2 9j — w\ — 3u72) + Fq cos(2 6 j — 4ru2) , (A2) 

where 

F° = fo + (e? + e|)/2 + el/4 + e?e |/5 + e\f& 

F° — eie 2 /io + e^/n + ei 62/12 
F* = e? el /17 

F2 = ei/45 + eiel/47 + ei/46 
Fl = eie2/49 + e?e2/5o + eiei /51 
Fo = 62/53 + 6162/54 + 62/55 
F 3 ~ e?e2/68 
Fi 1 = e^ei/eg 
Fl = e|/ 90 

Fl = e?e 2 /9i 

p2 _ 2 2 r 

-F 1 2 ~ ^1^2j92 

Fl = eiel /93 

Fl = el/ 94 . (A3) 

the quantities fi and fij are functions only of a = 
01/02 and their forms a re tabulated in Appendix B of 
iMurrav fc Dermottl d 19991 / 






